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A third-order Energy Stable Weighted Essentially Aon— Oscillatory (ESWENO) finite 
difference scheme developed by Yamaleev and Carpenter was proven to be stable in the 
energy norm for both continuous and discontinuous solutions of systems of linear hyper- 
bolic equations. Herein, a systematic approach is presented that enables “energy stable” 
modifications for existing WENO schemes of any order. The technique is demonstrated 
by developing a one-parameter family of fifth-order upwind-biased ESWENO schemes; 
ESWENO schemes up to eighth order are presented in the appendix. New weight func- 
tions are also developed that provide (1) formal consistency, (2) much faster convergence 
for smooth solutions with an arbitrary number of vanishing derivatives, and (3) improved 
resolution near strong discontinuities. 
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constant wave speed 

target value for w arising from stencil S( r ) 

matrix defining the discrete derivative operator 

artificial dissipation matrix operator 

skew-symmetric portion of derivative matrix 

symmetric skew-symmetric portion of derivative matrix 

energy-stable derivative operator 

discrete undivided difference operator 

continuous flux function (linear or nonlinear) 

discrete flux function built at position j + 1/2 

projection of continuous flux onto grid (linear or nonlinear) 

derivative of the flux function 

numerical flux function implicitly defining f(x) 

grid index 

symmetric positive definite matrix defining discrete norm 
skew-symmetric matrix defining dispersive portion of derivative 
symmetric matrix defining dissipative portion of derivative 
Stencil shifted one cell to the left, right 
solution 
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Weight function for stencil S( r ) 

Weight function for stencil <S( r ) , everywhere in the computational domain 
spatial coordinate 

solution dependent component of weight function arising from stencil S( r \ 

classical smoothness indicator for stencil r 

nonzero parameter used to prevent division by zero (ESWENO) 

grid spacing 

nonzero parameter used to prevent division by zero (WENO) 
diagonal matrix resulting from elemental decomposition of 5th-order 7 Z 
local component resulting from elemental decomposition of 1Z 
strictly positive local dissipation term 

fifth-order solution dependent component of weight function 
fifth-order solution dependent function proposed by Borges et al. 8 
parameter in one parameter family of 5th-order schemes 


I. Introduction 

A new, third-order weighted essentially nonoscillatory scheme (called Energy-Stable WENO) has recently 
been proposed and developed by the authors of the present paper. 1 In reference, 1 * we prove that the third- 
order ESWENO scheme is energy stable, that is, stable in an L 2 energy norm, for systems of linear hyperbolic 
equations with both continuous and discontinuous solutions. Stability is explicitly achieved (by construction) 
by requiring that the ESWENO scheme satisfies a nonlinear summation-by-parts (SBP) condition at each 
instant in time. Thus, L 2 strict stability is attained without the need for a total variation bounded (TVB) 
flux reconstruction or a large-time-step constraint, 23 and. 4 Herein, we generalize and extend the third-order 
ESWENO methodology 1 to an arbitrary order of accuracy. Similar to the third-order ESWENO scheme, 
the new families of higher order (up to eighth order) ESWENO schemes are provably stable in the energy 
norm and retain the underlying WENO characteristics of the background schemes. Numerical experiments 
demonstrate that the new family of ESWENO schemes provides the design order of accuracy for smooth 
problems and delivers stable essentially nonoscillatory solutions for problems with strong discontinuities. 

Another issue that is addressed in this paper is the consistency of the new class of ESWENO schemes. 
The consistency of any WENO-type scheme fully depends on a proper choice of the weight functions. On 
one hand, for smooth solutions the weights should provide a rapid convergence of the WENO scheme to 
the corresponding underlying linear scheme. On the other hand, the weights should effectively bias the 
stencil away from strong discontinuities. The high-order upwind-biased WENO schemes with conventional 
smoothness indicators that are presented in reference 2 are too dissipative for solving problems with a large 
amount of structure in the smooth part of the solution, such as direct numerical simulations of turbulence, 
or aeroacoustics, 5 . 6 Furthermore, as has been shown in references 7 and, 8 the classical weight functions of 
the fifth-order WENO scheme fail to provide the design order of convergence near smooth extrema, where 
the first derivative of the solution becomes equal to zero. New approaches are proposed in references 7 and 8 
to improve the error convergence near the critical points. Although these new weight functions recover the 
fifth order of convergence of the WENO scheme near smooth extrema, the problem persists if the first- and 
second-order derivatives vanish simultaneously. 8 An attempt to resolve this loss of accuracy is presented in 
reference. 8 This proposed resolution provides only a partial remedy for the problem; the same degeneration 
in the order of convergence occurs if at least the first three derivatives become equal to zero. To fully resolve 
this problem, we propose new weights to provide faster error convergence than those presented in reference, 8 
and impose some constraints on the weight parameters to guarantee that the WENO and ESWENO schemes 
are design-order accurate for sufficiently smooth solutions with an arbitrary number of vanishing derivatives. 

This paper is organized as follows. In section 2, energy estimates for the continuous and corresponding 
discrete wave equations are presented. In section 3, we present a one-parameter family of fifth-order WENO 
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schemes; one value of the parameter yields a central scheme that converges with sixth-order accuracy. 
In section 4, we present a systematic methodology for constructing ESWENO schemes of any order and 
demonstrate the methodology by transforming the family of WENO schemes presented in section 3, into a 
family of fifth-order ESWENO schemes. In section 5, we analyze the consistency of the new class of ESWENO 
schemes and we derive sufficient conditions for the weights functions that ensure that the ESWENO schemes 
are design-order accurate regardless of the number of vanishing derivatives in the solution. The tuning 
parameters in the weight functions are also optimized. In section 6, we present numerical experiments that 
corroborate our theoretical results. We summarize and draw conclusions in section 7. 

II. Energy Estimates 

Consider a linear, scalar wave equation 

I y + fy = °> / = an, t > °, ° < £ < t, 

u(0, x) = uq(x) 

where a is a constant and uq(x) is a bounded piecewise continuous function. Without loss of generality, 
assume that a > 0, and further assume that the problem is periodic on the interval 0 < x < 1. Applying the 
energy method to equation (1) leads to 

>IL = » m 

where || • \\l 2 is the continuous L 2 norm. Thus, the continuous problem defined in equation (1) is neutrally 
stable. 

We now develop using mimetic techniques (see 1 or 9 ) a class of discrete spatial operators that is neutrally 
stable or dissipative. The continuous target operator used for this development is the following singular 
perturbed wave equation: 


N 

S+S= Ef-ir-'igqefi*). /=«», 

n= 1 

w(0,x) = u 0 (x) , 


t > 0, 0 < x < 1, 


( 3 ) 


where n = fi(u) is a non-negative C°° function of u. As before, we assume that equation (3) is subject to 
periodic boundary conditions. Our goal is to match each spatial term in equation (3) with an equivalent 
discrete term that maintains neutral stability (or dissipates) of the discrete energy norm. 

We begin by showing that the terms on the right-hand side of equation (3) are dissipative thereby ensuring 
stability. Multiplying equation (3) by u and integrating it over the entire domain yields 


Id 
2 dt 


IMlL 


2 aU 
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= £ 

0 n= 1 ' 


/ \77 — I V ^ U, 

(— 1) u— — (i— — dx . 


dx 


dx r ' 


( 4 ) 


Integrating each term on the right-hand side by parts and accounting for periodic boundary conditions yields 
the following energy estimate: 


-I 

dr 


N 


= ~ 2 £ / 


a 


n — 1 


V dx r 


dx < 0 


( 5 ) 


All the perturbation terms included in equation (3) provide dissipation of energy. 

Turning now to the discrete case, we define a uniform grid Xj = j Ax, j = 0, J, with Ax = 1/J. 
On this grid, we define a flux f = au and its derivative f x = au x , where u = [u(xo, t), . . . , u(xj, t)) T 
and tO = [u x (xo, t), . . . , u x (xj,t)] T are projections of the continuous solution and its derivative onto the 
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computational grid. Next, we define a pth-order approximation for the first-order derivative term in equation 
(1) as 

OP 

— = Di + 0( Ax p ) . (6) 

ox 

Placing a mild restriction on the generality of the derivative operator (see 1 or 9 ), the matrix D can be 
expressed in the following form: 


D = P~ 1 [Q + R } ■ Q + Q T = 0 

R = R t ■ v T f?v > 0 (7) 

P = P T ; v T Pv > 0 

for any real vector v ^ 0. By choosing the matrix R as a discrete analog of the dissipation operator in 
equation (3), we have 

N 

R = Y J D i n M D i n } T , (8) 

n— 1 

where A is a diagonal positive semidefinite matrix and D\ is the difference matrix 


D i 


0 

-1 1 

0 


By using the SBP operators [eqs. (6)— (8)], the semi-discrete counterpart of equation (1) becomes 

fu - N 

— +P~ 1 Qi=-J2P~ 1 Di n A[D 1 n ] T i, (9) 

n— 1 

where f = au, u = [uo(t),Ui(t ), . . . , uj(t), ] T is the discrete approximation of the solution u of equation (1), 
and Q and A are nonlinear matrices (i.e. , Q = Q( u) and A = A(u). To show that the above finite-difference 
scheme is stable, the energy method is used. Multiplying equation (9) with u 4 1 P yields 

1 / N 

--\\uf P + a u T Qu = - a Y^ ([I?i n ] T u) T A[£)i"] T u, (10) 

” n— 1 

where || • || p is the P norm (i.e., ||u||p = u J Pu). Adding equation (10) to its transpose yields 
, N 

^||u||p + au T (Q + Q t ) u = -2 ([Hi”] 7 u) T A[£>i n ] 2 u . (11) 

n— 1 

If we account for periodic boundary conditions and the skew-symmetry of Q , then the second term on the 
left-hand side vanishes, and the energy estimate becomes 

|]u||7> = -2a f; ([H 1 ”] T u) T A[H 1 n ] T u < 0 . (12) 

n— 1 

The right-hand side of equation (12) is nonpositive because the diagonal matrix A is positive semidefinite 
[v T Av > 0 for all real v of length ( J + 1)] and a > 0; thus the stability of the finite-difference scheme given 
by equation (9) is assured. This result can be summarized in the following theorem: 

Theorem 1. The approximation [eq. (9)] of the problem [eq. (1)] is stable if equations [(6)-(8)] hold. 
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Remark 1. Despite the fact that the initial boundary value problem [eq. (1)] is linear, the finite-difference 
scheme [eq. (9)] constructed for approximation of equation (1) is nonlinear, because the matrices Q and A 
(and in principle P) are assumed to depend on the discrete solution u. 

Remark 2. The only constraints that are imposed on the matrix Q and the diagonal matrix A are the skew- 
symmetry of the former and positive semidefiniteness of the later. No other assumptions have been made 
about a specific form of the matrices Q and A to guarantee the stability of the finite-difference scheme [eq. 

( 9 )]- 

Remark 3. The discrete operators that are defined by equations [(6)— (8)] are similar in form to those that 
are used for conventional SBP operators (See refs. 9,10 ). What is new, however, is the fact that the matrices 
Q and A depend on u. 
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Figure 1. Extended six-point stencil S' 6 , and corresponding candidate stencils Sll , Sl , Sr, and Srr for one-parameter 
family of fifth-order WENO schemes. 


Next, a new one-parameter family of fifth-order WENO schemes is developed, and then is used as the 
starting point in the development of a family of “energy-stable” WENO schemes. Great care is exercised in 
the developing the WENO schemes to ensure that design-order accuracy is achieved in the vicinity of smooth 
extrema. 


III. Fifth- and Sixth-order WENO Schemes 


Any conventional high-order WENO finite-difference scheme for the scalar one-dimensional wave equation 
(1) can be written in the following semidiscrete form: 


du 3 , fjj± l±± = n 

dt Ax 


(13) 


For the fifth-order WENO scheme that is presented in reference, 2 the numerical flux /- + 1 is computed 
as a convex combination of three third-order fluxes defined on the following three-point stencils: Sll = 
{xj- 2 ,Xj-i,Xj}, Sl = {xj-i,Xj,Xj + i}, and Sr = {xj,Xj+i,Xj+ 2 }- (See figure 1.) Note that this set of 
stencils is not symmetric with respect to the ( j + |) point; thus, the fifth-order WENO scheme is biased 
in the upwind direction. A central WENO scheme can be constructed from the conventional fifth-order 
WENO scheme by including an additional downwind candidate stencil Srr = {xj+i, Xj+ 2 , Xj+ 3 }, so that 
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the collection of all four stencils is symmetric with respect to point (j + 4). a The WENO flux, constructed 
in this manner, is given by 



„LL 


cLL 


~ W j+l/2Jj+l/2 + W j+l/2Jj+l/2 


+ W 


R fR 
j+l/2Jj+l/2 


, HH fHH 
+ W j+l/2Jj+l/2 


where /j+j/ 2 > r = {LL, L , R, RR } are third-order fluxes defined on these four stencils: 


/ f ( Uj - 2 ) \ 


( f LL {u j+ 1 / 2 ) \ 


^2 -7 11 0 ^ 


f(uj- 1 ) 

f L ( u j+ 1 / 2 ) 

1 

-15 2 


f{Uj) 

f R (u j+ 1 / 2 ) 

6 

2 5-1 


f{Uj+ 1 ) 

V f RR {u j+ 1 / 2 ) ) 


^ 0 11 -7 2 ) 


f(uj+ 2 ) 


V /( w f+s) / 


(14) 


(15) 


and w LL , w L , w R , and are weight functions that are assigned to the four stencils Sll, Sl, Sr, and Srr, 
respectively. The terms w LL , w L , w R , and w RR in equation (14) are nonlinear weight functions. These 
have both preferred values that are derived from an underlying linear scheme as well as solution-dependent 
components. The preferred values are given by 


d LL = ± j-p; d L =^~ 3p-, d R = ± + 3p ; d RR = p , (16) 

where p is a parameter. The convergence rate of the scheme [eqs. (13) — (16)] with the preferred values 

u ,( r ) 

= d ^ and r = {LL, L, R, RR} is greater than or equal to 5 for all values of the parameter ip in equation 
(16). For the specific value <p c = 75^, the fifth-order term vanishes and the convergence rate increases to 6. 
The classical fifth-order upwind-biased WENO scheme of Jiang and Shu is obtained for p = 0. 

The weight functions w LL , w L , w R , and w RR needed in equation (14) are given by 



E< 


(17) 


where 

a r = d M ^1 + , r = {LL, L, R, RR}, (18) 

and e < 0( Ax 2 ). The functions are the classical smoothness indicators. [See equation (59) for the 
general expression.] For fifth-order schemes, they are given by 


(3 LL = 

+Uj) 2 

+ 

\ (Uj - 2 - 4Uj_i + 3 Uj) 2 

P L = 

12 ( U f— 1 _ 2 11 j + Uj+l) 

+ 

\ (Uj-! -Uj+l) 2 

1 3 R = 

12 (' u j ~ 2uj + 1 + 2 ) 

+ 

\ (3 Uj - 4u i+ i + u j+2 ) 2 

P R R = 

12 ( U j + 1 — 2Uj_|_2 + Uj+ 3 ) 

+ 

\ (-5tij+i + Suj+z - 3 u j+3 ) 


The T p also varies with discretization order; the expression for 75 is given by 


T5 


i~fj- 2 + 5/j-i - 10 fj + lOfj+i - 5 f j+2 + f j+3 f , for ip ± 0 
(fj - 2 - 4/j-i + 6 fj - Af j+1 + f j+2 f , for p = 0 


(20) 


Expressions for the fourth-, seventh-, and eighth-order WENO schemes appear in the appendix. 

a The above approach to constructing central WENO schemes is proposed in reference. 6 In general, central high-order WENO 
schemes that are built in this manner, are unstable when unresolved features or strong discontinuities exist in the computational 
domain. The generalized energy-stable methodology presented in the next section guarantees the stability of the even order 
approximations while maintaining their nonoscillatory properties. 
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The WENO mechanics expressed in equations (17)-(20) represents a significant departure from the 
mechanics used in the original algorithm by Jiang and Shu. 2 For example, equation (18) replaces the 
expression 


d (r) 


(21) 


that is used in the weight functions given in reference. 2 Furthermore, e < 0( Ax 2 ) replaces e = O(10 -6 ) in. 2 
Equations (17)-(20) more closely resemble the weight and smoothness indicators proposed by Borges et al. 
in 7 8 There are differences, however, in how the parameters t p and e are chosen; differences motivated by the 
following observations. 

Borges et al. 8 proposes the following function t 5 : 


T5 = \P LL - P R \ 


(22) 


and a fixed value for the parameter e. Although the above weights and smoothness indicators [eqs. (17)— (19)] 
with the value of T 5 given by equation (22) significantly outperform the conventional weights of Jiang and 
Shu 2 for both continuous and discontinuous solutions, three remaining shortcomings include: 

1. The value of given by equation (22) is not smooth because an absolute value function is used, which 
may reduce accuracy at points where (3 LL — (3 R changes sign. 

2. This approach currently does not generalize to WENO schemes with design orders other than 5. For 
example, neither tq = \/3 LL — (3 R \ nor tq = \(3 LL — fi RR \ provides the design order of accuracy for the 
central sixth-order WENO scheme. 

3. Near critical points where /', f", and f" vanish simultaneously, the modified weights [eqs. (17)- 
(19), and (22)] fail to provide the design order of convergence for the fifth-order WENO scheme if no 
constraints are imposed on the parameter e other than e > 0. An example that demonstrates this 
property is presented in Section 6.1. 

Selecting 75 in equation (20) to be a quadratic function of the fifth-degree undivided difference defined on 
the entire six-point stencil, circumvents the first, and second drawbacks encountered when using the scheme 
suggested in reference 8 . Indeed, 75 is a C°° function in its arguments, and can readily be generalized to 
WENO schemes of order p by choosing t p to be the pth-degree undivided difference that is defined on the 
entire (p + l)-point stencil. In contrast to the T 5 found in reference, which is of order 0(Ax 5 ) for smooth 
solutions, the proposed function 75 as given by equation (20) is of order 0( Ax 8 ) and 0(Aa; 10 ) for ip = 0 
and <p 0, respectively; thus much faster convergence of the fifth- and sixth-order WENO schemes to the 
corresponding underlying linear schemes is achieved. b 

A remedy for the third shortcoming requires an additional modification to the approach proposed in 
reference . 8 Both the fifth- and sixth-order WENO schemes with the new weights that are given by equations 
(17-20) are design-order accurate for smooth solutions, including points at which the first- and second-order 
derivatives of the solution vanish simultaneously. However, if all derivatives up to fourth are equal to zero, 
then the fifth- and sixth-order WENO schemes locally become only third-order accurate. To fully resolve 
this issue, the constraint e < 0(A:r 2 ) is imposed herein. 

Equations (13)-(20) describe a new family of fifth-order WENO schemes. The primary difference between 
existing WENO schemes and this new family, is in the stencil biasing mechanics described by equations (17)- 
(20). Detailed theoretical justifications for the choice of t p and e used in equations (17)-(20) are presented 
in sections 5.2, 5.3. There, and again in the results section, it is shown that these parameters provide fast 
convergence of the new WENO and ESWENO schemes to the corresponding underlying linear schemes for 
smooth solutions, and deliver improved shock-capturing capabilities near unresolved features. 

b The conventional fifth-order WENO scheme that corresponds to -p = 0, does not include the downwind stencil Srr; 
therefore, the fourth-degree undivided difference built on the available five-point stencil, is included in equation (20). This 
approach to handle the narrow stencil encountered with if = 0, also generalizes to other orders of accuracy. 
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IV. A General Approach to Constructing High-Order Energy-Stable Schemes 

Algorithm (1) transforms an existing WENO scheme into an ESWENO scheme that is characterized 
by the following four properties: 1) a bounded energy estimate for arbitrary nonsmooth initial data, 2) 
conservation, 3) design order accuracy for sufficiently smooth data, and 4) discontinuity (shock) capturing 
capabilities that are similar to those of the base WENO scheme. 

Algorithm 1 ( Transformation from WENO to ESWENO). 


1 : 

2 : 

3: 

3a: 


3b: 


3c: 

4 • 


Express the base WENO derivative operator as matrix D. 

Decompose D into symmetric and skew- symmetric components as D = D s k + D s ym 
Add an artificial dissipation operator D a d such that the modified symmetric matrix 
D sym + D a d is positive semidefinite. 

Form the decomposition D sym = ^®_ 0 [-Di*]Aj[.Di*] > 

2s + 1 is the bandwidth of the matrix D. 

The existence of this decomposition is established in the appendix. 

Modify the diagonal terms of A,; such that they are smoothly positive. 

One such approach is 

(A* )j,jl 


(A i)j,j ~ 2 


(Ad 2 


jj 


6 ? 


where Si, 0 < i < s are small positive constants that may depend on Ax. 
Form the symmetric, positive semidefinite matrix 

D ad = P^ELolDilMDif 

Form the energy-stable operator 

D = D + D a d 


(23) 


(24) 

(25) 


By construction, algorithm (1) is applicable to 1-D periodic schemes of any order and produces a modified 
scheme that automatically satisfies the first and second properties: bounded energy estimate and conservation 
(See reference 1 for details.) Likewise, as shown in the next section, the additional terms added in the 
ESWENO formulation do not degrade the formal accuracy of the original WENO discretization (property 
three). However, the degree to which the two formulations differ for unresolved data is not clear. Thus, the 
properties of the new ESWENO scheme must be tested to ensure it retains the desirable properties of the 
original formulation. 

ESWENO Schemes of Fifth- and Sixth-order 

We now apply algorithm (1) to the one-parameter family of fifth-order WENO schemes [eqs. (17)— (20)] for 
the scalar 1-D wave equation (1) with a > 0. Combining equation (15) with equation (14) and substituting 
the resulting WENO flux / J+ 1 into equation (13) produces a stencil for the jth grid point of the form 


D 


5 


/ 0 

2 ”JA/2 

- 7 ™j+l/2 

11 ~t+l/2 

0 

0 

0 

\ 

0 

0 

- 1W j+ 1/2 

5 “t + l/2 

2 ™j + l/2 

0 
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0 

0 

0 

H + l/2 

5 ”j+l/2 

-1 “J + l/2 

0 


0 

0 

0 

0 

1 1 luRR 
llw j- |-l/2 

7w j+ 1/2 

2w j + 1/2 



^-1/2 


0 

0 

0 

0 


0 

1W i- 1/2 

~ 5w j- 1/2 

~ 2w j‘-l/2 

0 

0 

0 


0 

0 

- 2 ”f-l/2 

“ 5 ”f-l/2 

1 “f-l/2 

0 

0 


V 0 

0 

0 


7 3- 1/2 

- 2 “f- B l/2 

0 

) 


\ 


) 


( 26 ) 


with k on the interval j — 3 < k < j + 3. An explicit expression for the differentiation matrix D 5 follows 
immediately from equation (26). 
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The derivative matrix D 5 is now decomposed into symmetric and skew-symmetric parts as 


n-5 _ t ) 5 

x ^ skew 


D 


sym • 


As with the third-order case (ref. 1 ), the skew-symmetric component of D 5 and the norm, take the forms 


D° 


= P~ l Qb ; Q 5 + Qi = 0 ; P = Axl. 


while the matrix D 5 sym is expressed as 

D%m = P- 1 ( D\ A® [Dlf + D\ k\ [Dl\ T + D\ A\ [D\\ T + £>? Aq [£>?]’ 
The matrices A^ are diagonal with expressions for the jth element dehned by 


(27) 


( A 3 )j,j ~ 


<+<-<+< 


(28) 


(A 2 ) 


hi 


12 


+<+3/2 - Aw f+ 5/2 + <+3/2 
-<+ 1/2 + 4 <-l/2 - <+1/2 . 


(A?) 


3,3 




‘XnnLL 
6W j+ 1/2 

- 5 «2 

+2 <+5/2 

1 


~ >rW j+ 1/2 

-<+3/2 


12 

+<-1/2 

_ <+l/2 




_ -2« 2 +5« 2 

_3 v) RR 
6W j+ 1/2 





1 

2 


~ W j-l/2 ~ W j- 1/2 

+<+ l / 2 + <+ 1/2 


- <_ 1/2 - « 2 
+ <+ 1/2 +<+<2 . 


= 0 


(29) 


(30) 


(31) 


Because = 1> is always equal to zero, and is not included in D\ d . The sign of the diagonal terms 

( A '' 5 )+j> * = 2, 3, could be either positive or negative; thus, the conventional fifth- and sixth-order WENO 

schemes may become locally unstable. If we define (A f ) • ■, i = 1,2, 3, to be smoothly positive 


(Af< = . 

then the additional artificial dissipation operator becomes D b ad = P _1 ^®_i[.Di*](Af )[£>!*] , and the re- 
sulting energy-stable scheme is obtained by adding the additional dissipation term to the original WENO 
scheme. That is, 

D 5 = D 5 + D 5 ad . (32) 

By construction, D 5 satisfies all of the conditions of Theorem 1, thereby providing stability of the fifth- and 
sixth-order ESWENO schemes that are defined by equation (32). 


V. Consistency Analysis 

V.A. Necessary and Sufficient Conditions for Consistency of WENO Schemes 

(r) 

We now derive the necessary and sufficient conditions for the weight functions w; ’ for a family of pth- 
order WENO schemes to attain the design order of accuracy. Similar conditions have been obtained for the 
conventional fifth-order WENO scheme in reference. 7 With the same approach discussed in section 3 for 
p = 5, a one-parameter family of ptlr-order WENO fluxes can be constructed by using a convex combination 
of (s + 1) fluxes /< of sth order as 

fj±i/2 = <+i / 2 /+H /2> (33) 
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with 


fj± 1/2 = h (Xj± 1/2 ) + XM r)A ‘ T * + °( Aa;P+1 ) 


(34) 


l=S 

where h(a;) is the numerical flux function that is implicitly defined as 

x+ ^ 


f(x) = 


1 

Ax 


Hv ) dr] 


(35) 


and c\ ' are constants that do not depend on Ax. 


The corresponding pth-order WENO operator that approximates the first-order spatial derivative is given 


by 


mi = 


fj+ 1/2 - fj~ 1/2 


?(■ 


(r) /»(r) 

w j+i/ifj+ 1/2 - w 


(r 

j-l/Uj- 


f (r) i 

h- 1 / 2 ; 


(36) 


J Ax Ax 

where [■] ,■ is a jth component of a vector, the index r in equation (36) sweeps over all (s + 1) stencils, and 
is a nonlinear weight function that is assigned to the corresponding s-point stencil S r . For sufficiently smooth 
solutions, the weights 't±l T> approach their preferred values d^ r \ so that the WENO operator converges to 
the target linear operator £) Tar « et as 


fT-get _ .Target £ (d« f^ /2 - d« L/2 ) 

j-^Targetfj = 0+1/2 £-1/2 _ r V ^ 7 ■> ' J 


Ax 


Ax 


(37) 


where d 7 ' -* is a one-parameter family of constants chosen to ensure pth-order convergence of the target 
operator to the exact value of the first-order derivative at Xj. That is, 


[Z) Target f = 


df 

dx 


0{ Ax p ). 


(38) 


The coefficients d 64 for one-parameter families of linear target schemes up to eighth order are given in the 
appendix. All target linear schemes in the family are (2s — l)th-order accurate, except one central scheme, 
which is ( 2 s)tli-order accurate. 

Subtracting equation (37) from equation (36) and using equation (34), we have 


E[(™'+ ) i/2 - dW )/i+\ /2 -(™)- 1/2 - d<r) )/i- 1/2 . 


Ax 


[D p {\ ■ - [£)Target f j . = 

= Mi/2 - d (r) )h (r) (x j+ 1 /2 ) - {uA l/2 - d^)h (r) (xj-i/2) 

r L 

+ E E Ax 1 - 1 ^ [(ii /2 - d (r) ) - (ii /2 - d (r) )l + 0( Ax p ) 


(39) 


From equation (39) it immediately follows that to retain pth-order accuracy, the weights of the WENO 
operator D p should satisfy the following necessary and sufficient conditions: 


E 


£’4 _ d(r) 

£± 1/2 “ 


E4 r) (■ 


= 0( Ax p+1 ) 

- dW) - (w { l l } 1/2 - dW)l = 0( Ax p ~ s+1 ) 


j+l/2 “ > ^3- 1/2 'j ' ( 40) 

E c P ) \( w j r +i/2 - d [r) ) ~ (i-1/2 “ d {r) ) 1 = °( A d 

r J 

Here, we use the following properties of h(x) and d ^ : f'(x) = ^ anc j ^^(r) _ 4 
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By construction, the weights [eq. (17)] are normalized such that J2 r ■ = 1; thus, the first constraint in 

equation (40) is satisfied identically. To simplify the analysis, especially for f(x) with an arbitrary number of 
vanishing derivatives, we hereafter use the following sufficient condition on w' for the conventional WENO 
scheme to attain pth-order accuracy: 

w ( r )-d (r) =0( Ax p ~ s+1 ) . (41) 

This constraint is a direct consequence of the necessary and sufficient conditions [eq. (40)]. 

V.B. Sufficient Conditions for Consistency of ESWENO Schemes 

In this section, we show that the conditions [eq. (41)] and the constraints on <5, in equation (23): 

Si = 0(Ax p ~ 2i+1 ) , 1 <i< a, (42) 

guarantee that the energy-stable modifications of the conventional pth-order WENO scheme (see section 
4) preserve the design order of the original scheme. As follows from equation (25), the ESWENO operator 
consists of two terms: one is the original WENO operator and the other is the additional artificial dissipation 
operator D a d that is given by equation (24). As shown in the previous section, the WENO operator is pth- 
order accurate if equation (41) holds. Hence, we need only show that D a di = 0( Ax p ). 

Let us prove this conjecture for the one-parameter family of fifth-order ESWENO schemes that is pre- 
sented in Section 4.1. Note that the term P~ 1 Df A 0 [£>f] T in equation (27) is identically equal to zero because 
of the normalization y) w- r ^ = 1; it need not be included in D a d . Thus, the additional artificial dissipation 
term is given by 

D ad f = P- 1 DiA 1 [Zl 1 ] T f + p- l DlK 2 [Dl] T i + p- 1 P?A 3 [D3] T f (43) 

A* = - \J A, 2 + 5f + Xi , * = 1,2,3 (44) 

where A i is a jth diagonal element of the matrix A,. For simplicity, we have omitted the superscript 5 and 
the subscript (j,j) in this section. 

First, we evaluate the terms Ai, A2, and A3 that are defined by equations (28-30). To simplify the 
derivation, the sufficient conditions [eq. (41)] for the one-parameter family of fifth-order ESWENO schemes 
is rewritten in the following form: 

w : (r) - d (r) = (v - °( Ax 3 ) + 0( Ax 4 ) . (45) 

In equation (45), the stencil width s of each reconstruction polynomial is equal to 3, and the order p of this 
family of schemes is equal to 5 if ip ^ A, or 6 if p = As follows from equation (45), the stiffer constraint 
on the weights should be imposed to obtain sixth-order accuracy. 

Replacing w[ T> in A2 with the corresponding preferred values SD and using equations (16) and (29), for 
any value of the parameter p yields 

^ [d LL - 4d Li + d L - d R + 4 d RR - d RR ] = 0 . (46) 

Subtracting equation (46) from equation (29) and taking into account equation (45) yields 

A2 = j 2 (™f+ 3/2 - dLL ) - A ( w f+, 5/2 - dLL ) + K L +3/2 - d L ) 

~( w f+i/2 - dR ) + 4 ( w f-i/2 - d RR ) ~ (wf+1/2 - d RR ) (47) 

= °( Ax 3 ) + Q{ Ax A ) . 
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By comparing equations (29) and (30), one can see that Ai is at least one order higher than that of A 2 
because of an additional cancellation that occurs within each group of terms that is associated with the same 
stencil. For example, expanding all of the terms w RL in equation (30) about Xj yields 


3 Wj+X/2 - 5u b+3/2 + 2 «b+ 5 /2 = 


dw LL 

dx 


Ax + 0( Ax 2 ) 

Xj 


(48) 


which gives an extra factor 0( Ax) compared with the corresponding terms w LL in equation (30): w RR 3 , 2 — 
4 w f+ 5/2 = ^(l)- The same conclusion can be drawn for the other groups of terms that are associated with 
the L, R, and RR stencils in equation (30), which leads to 


Ai 



0{Ax a ) + 0{Ax 5 ) . 


(49) 


Applying the same procedure to A 3 yields 


A3 


(■ w LL 


3+ 5/2 


- d LL ) - {w RR 


1 ( d LL _ d RRj + 1 

5 (lo _ 2( p) + W ~ to) °( Ax3 ) + 0(Ax A ) 


3+ 1/2 


d RR ) 


(50) 


To evaluate each term in equation (43), we consider three cases: 1) |Aj| » <5, > 0, 2) \ = 0{8i), and 3) 
| A,; | <C Si, i = 1,2, 3. If we assume that |A, | Si > 0, then equation (44) can be expanded as follows: 


which yields 


\ _ I'M + M 

‘ ” 2 + 4|Aj| ’ 


f Ai , if A i > 0 

A* = < Si , if Aj = 0 , 

{ if < 0 


(51) 


(52) 


where the higher order terms have been omitted. Because equation (52) has been derived under the as- 

s 2 

sumption that |Aij <5, > 0, we can immediately conclude that | Ai | Si ^ jpH-. Therefore, we only need 

consider that 


Ai — Ai 


(53) 


which provides the lowest order of convergence for the term D\Ai[D\\ T f. If we substitute equation (53) in 
equation (43) and use equations (47, 49, 50), then the additional ESWENO dissipation term becomes 


D ad f = [{'p- A 5 )0{Ax 3 ) + 0{Ax i )}D 1 [D 1 } T i 

+ [(^ “ to) °( Ax 2 ) + 0(Aa; 3 )] Dl[D\] T i 


+ 


To- 2 ^ 

6Ax 


+ fa - to) 0(Ax 2 ) + Q(Ax 3 ) D\[Dl \ T f 


If we take into account that D\[D \] T f = 0{ Ax 21 ), then D a di can be recast as 

Dad f =(e-^) °( A x 5 ) + 0( Ax 6 ) . 


(54) 


(55) 


From equation (55) we see, that the additional dissipation term is at least fifth-order accurate for all values 
of the parameter ip. The fifth-order term vanishes for ip = A.; for this value the order increase by one to 
sixth-order. 

The second case can be considered in a similar manner. Substituting A * = 0(Si), i = 1,2,3, in equation 
(43) yields 


Dadi 


°gLLD 1 [D l ] T f + °g£iDl[Dl) T f + ^D\[Dl} T i 
0(5\Ax) + 0(S 2 Ax 3 ) + 0(S 3 Ax 5 ) . 


(56) 
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To guarantee that the ESWENO dissipation term is pth-order accurate, the following constraints should be 
imposed on d t : 

51 = OiAxP- 1 ) 

5 2 = 0(AxP~ 3 ) (57) 

53 = 0(AxP~ 5 ) , 

where p is equal to 5 for ip ^ ^ 0 r 6 for ip = Note that the constraints [eq. (57)] are fully consistent 
with those that are given by equation (42). The third case, |Ai| <C Si, is similar to the second one and results 
in the same constraints [eq. (57)] on 5,; therefore the third case is not presented here. 

Remark 4. Note that Si, i = 1,2,3, are user-defined parameters; therefore, the conditions [eq. (42)] can 
always be met. 

Remark 5. Although only fifth- and sixth-order ESWENO schemes have been analyzed in this section, the 
same procedure is directly applicable to the additional ESWENO schemes presented in the appendix. Thus, 
we can conclude that if equations (41) and (42) hold, then all the ESWENO schemes that are considered in 
this paper are design-order accurate. 


V.C. Consistency of the ESWENO scheme with New Weights 


The new weight functions for the one-parameter family of pth-order WENO and ESWENO schemes are 
given by 


Qlr- 


w y ' = 




a , 


= d (r) ( 1 + 


+ /3( r 


(58) 


where ft 1 ' 1 ' 1 are the classical smoothness indicators: 


/3 (r) 


1=1 



{ d l q r (x) 
\ d l x 


2 

dx , 


(59) 


q r (x) is an (s — l)th-degree reconstruction polynomial that is defined on a stencil S r , and e is a small positive 
parameter that can depend on Ax. In equation (58), r p is defined by 


t p = (V < Xj- S+ 1 , . . . , x j+s >) , for <p y 2 0 (60) 

t p = (V < Xj - S+ 1 , . . . , xj +s - 1 >) 2 , for <p = 0 (61) 

where V < Xj- s+ \, . . . , Xj +S > is the pth-degree, undivided difference. Note that for the original WENO 
schemes of Jiang and Shu , 2 which correspond to p = 0, the entire stencil includes only (2s — 1) points; 
therefore, the highest degree of undivided differences that can be constructed on this stencil is 2 s — 2 , rather 
than 2s — 1, which is used for the other schemes in this one-parameter family. In particular, w^ r \ /3^ r \ and 
T 5 for the fifth- and sixth-order WENO and ESWENO schemes are given by equations (17)-(20). Another 
scheme that requires special consideration is the central ESWENO scheme, which is obtained by setting 
p = ipc, so that it is one order higher than that of the other schemes in the family. 

First, we present a truncation error analysis for the entire one-parameter family of pth-order ESWENO 
schemes, except for the two schemes that correspond to ip = 0 and p = p c . We demonstrate that the 
new weights that are defined by equations (58)-(60) satisfy the sufficient condition [eq. (41)] for smooth 
solutions with any number of vanishing derivatives if the following constraint is imposed on the parameter 
e in equation (58): 

e > 0(Ax p+s_1 ) > 0, (62) 

where p is a design order of the scheme and (s — 1 ) is a degree of the corresponding reconstruction polynomials. 

If we assume that all of the required derivatives are continuous and use the properties of the Newton 
divided differences, then the Taylor series expansions of (3^ and t p (<p ^ 0) at Xj are given by 


/? (r) = f 2 Ax 2 + 0(Ax 3 ) , 


( 63 ) 
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Tp = (/ (p) ) 2 Ax 2p + 0(Ax 2p+1 ) , (64) 

where f and / ^ are the first- and pth-order derivatives of / at Xj. For example, for the family of fifth-order 
WENO schemes, these expansions are 

(3 LL = f 2 Ax 2 + (if/" 2 - Iff") Ax 4 + (-f f'f" + Iff") Ax 5 + 0{ Ax 6 ) 

/3 l = f 2 Ax 2 + (g/" 2 + Iff") Ax 4 + 0( Ax 6 ) 

P R = f 2 Ax 2 + }^f" 2 - Iff") Ax 4 + (f - Iff") Ax 5 + 0(Ax e ) 

pRR = f 2 Ax 2 + ( y §/" 2 - fff") Ax 4 + (fff - 5/7"") A;r 5 + 0( Ax 6 ) 

T5 = Aa- 10 + 0( Ax 11 ) , for <p f 0 . (66) 

We first consider a case with f(xj) f 0. Substituting equations (63) and (64) in equation (58) and 
accounting for equation (62) yields 


= 0 ( Ax 2p ~ 2 ) (1 - 0(Aa; p+s_3 )) , (67) 

which leads to 

«,(*■)= dW + 0( Ax 2p ~ l ) . (68) 

Note that the order of convergence of uf to its preferred value is 2p — 1 rather than 2 p — 2. The main 
reason for such “superconvergence” is the additional cancellation that occurs because the leading truncation 
error terms of all of the smoothness indicators are identical to each other if f f 0, as can be seen in 
equation (63). Equations (67) and (68) are valid only for p > 3 and are not applicable to the third-order 
WENO and ESWENO schemes that correspond to p = 0. The detailed analysis of these third-order schemes 
(ip = 0) is presented in reference. 1 By comparing equation (68) with equation (41), we can immediately 
conclude that the new weights satisfy the sufficient condition [eq. (41)], ensuring that both the WENO and 
ESWENO schemes are design-order accurate. Furthermore, for schemes of order 3 (<p / 0) or higher, the 
weights converge to their preferred values at a rate that is significantly faster than that given by the sufficient 
condition (41). The result is a much faster rate of convergence for both the ESWENO and WENO schemes 
to the corresponding target linear schemes, even on coarse and moderate grids. 

The next issue that is addressed is the convergence of the ESWENO schemes near the critical points at 
which the first-order and higher order derivatives of the flux approach zero. Let x c be a critical point at 
which the flux function is sufficiently smooth and at which its derivatives of order up to n v d th are equal to 
zero. That is, 

f(x c ) = ■■■ = f M (x c ) = 0 , f nvd+1 \x c ) f 0 . 

In contrast to the previous case for which f f 0, the leading truncation error terms of the smoothness 
indicators at the critical point are not equal to each other; thus, no additional cancellation occurs. For any 
number of vanishing derivatives, the following inequalities always hold: 

e > 0(Acc p+s_1 ) > 0( Ax 2p ) > r p , 


which yields 


By using equation (69), the weights 


e + 

can be recast as 

1 + T,7W t ' 


« 1 . 


d (r) + o 


T P 
6 + 



(69) 


(70) 


14 of 30 


American Institute of Aeronautics and Astronautics 



where we use d ^ = 1. For any number of vanishing derivatives, we have 


: + f3( r 


<^< 


0(Ax 2p ) 

0(Ax p+s ~ 1 ) 


= 0( Ax p ~ s+1 ) . 


( 71 ) 


From equation (71) we see that w ^ converges to d ^ at the rate of 0 (A.t p_s+1 ) or higher and satisfies the 
sufficient condition (41). 

As mentioned at the beginning of this section, the two schemes that correspond to <p = 0 and <p = <p c 
require special consideration. Using the same procedure that is outlined above, we can easily show that if 
the parameter e satisfies the constraints 


e > 0(Aa: 3s_4 ) , for p = 0 
e > 0(Ax 3s ~ 3 ) , for ip = p c , 

then the corresponding ESWENO schemes are design-order accurate regardless of the number of vanishing 
derivatives of the solution. The constraints [eq. (72)] are derived using the following relations between the 
order of the scheme and the degree of the reconstruction polynomials: 


p = 2s — 1, for (p = 0 
p = 2s, for ip = ip c . 


Also, note that if no constraint is imposed on e except that it must be strictly positive, then the order 
of convergence of the pth-order WENO and ESWENO schemes with the new weights may deteriorate from 
p to s. Indeed, for a sufficiently large number of vanishing derivatives n v d, (3^ and r p may become of the 
same order. For example, for the family of fifth-order schemes (^ / 0), this type of degeneration occurs at 
n VC L = 4. If /' = f" = /"' = f"" = 0, f( 5 ' 1 ^ 0, and e < 0(Aa; 10 ), which does not satisfy the condition given 
in equation (62), then equations (65) and (66) lead to 


T5 

e + (3 M 


0(Aa: 10 ) 

0(Ax w ) + 0(Ax 10 ) 


0 ( 1 ) . 


(73) 


As a result, the weights are of order 0(1) and the fifth-order WENO and ESWENO schemes locally degen- 
erate to third order. 

Remark 6. In reference, 8 the following modification of the weight functions given in equations (17) — (19), 
and (22) is proposed to recover the fifth-order rate of convergence if n v d < 2: 


Jr) — 


l 


= d (r 


i + 


p 5 


-P(r) 


(74) 


where m = 2. However, the modified weights [eq. (74)] still experience the same degeneration in accuracy 
near critical points for any choice of the parameter m if n v d < 6. We demonstrate this with the following 
“thought” experiment. First, build a polynomial of degree 7 (or higher) such that at least 6 derivatives vanish 
at a single point. Next, connect this polynomial to a constant polynomial at the point where the derivatives 
of the first polymonial vanish. The solution to equation (1) is then a six times continuously differentiable 
function everywhere on the combined piecewise polynomials. Note that the underlying linear scheme is 
fifth-order accurate in this case. If we assume that the LL stencil is located completely in the constant part 
of the solution (i.e. fj -2 = fj-i = fj ) while the other two stencils include points from the polynomial part 
of the solution, then we have (3 LL = 0, f3 L ^ 0, and (3 R ^ 0. As a result, 75 = \(3 LL — P R \ = f3 R . If no 
constraint is imposed on e, then we can always choose e such that j3 R > £ -» 0 on any given grid, which 
yields 

This leads to w R = 0(1) and, consequently, to a loss of the design order of accuracy. Note that this problem 
persists regardless of the choice of m in equation (74). 
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Remark 7. The parameter e is user-defined, and therefore, the sufficient conditions [eqs. (62) and (72)] can 
always be satisfied. 

Remark 8. Equations (62) and (72) do not provide sharp estimates for the parameter e, which can be weak- 
ened if additional information regarding the number of vanishing derivatives is available a priori. Note, 
however, that the constraints [eqs. (62) and (72)] guarantee the design order of convergence of the cor- 
responding WENO and ESWENO schemes for smooth solutions with an arbitrary number of vanishing 
derivatives. 

The last issue that we discuss in this section concerns discontinuous and unresolved solutions. To suc- 
cessfully emulate the ENO strategy, a stencil for which the solution is discontinuous is eliminated from the 
approximation by effectively nullifying the corresponding weight that is associated with this stencil. Let a 
discontinuity be located inside a stencil S r , while the solution is smooth in all other stencils. We can easily 
verify that t p = 0(1), because t p involves all points of the entire stencil including those that contain the 
discontinuity. Therefore, the weight can be evaluated as 


w 


M _ 


d (r) 


1 + 


o(i) + E 

l^r 


0 ( 1 ) 

e+0(l) 

0 ( 1 ) 

e+0( Ax' 2 ) 


0 ( 1 ) 


0 ( 1 ) 


0 ( 1 ) 

e+0( Ai 2 ) 


0(1) [e + 0( Ax 2 }) 


(75) 


Based on the above equation, to reduce the influence of e on the solution near the discontinuity, the parameter 
e should satisfy the following constraint: 

e < 0(Ax 2 ) . (76) 

Indeed, if equation (76) is met, then w ^ is of the order 0(Ax 2 ) and has the same order of magnitude as it 
would have if e = 0. Hence, the parameter e must be bounded not only from below by equations (62) and 
(72), but also from above by equation (76). 

Another consideration that can help us optimally select the parameter e is the manner in which the 
WENO and ESWENO schemes handle small-amplitude oscillations. Consider a solution that contains small- 
amplitude spurious oscillations 

f = fs + Sf , (77) 

where f s is a smooth component and 5f is a nonsmooth, high-frequency component of the solution. By 
substituting equation (77) in equation (59) and assuming that the contribution of the smooth component is 
negligibly small, we have 

/3 (r) = 0(5 f 2 ) . 

As follows from equation (58), if e 7§> 0(5 / 2 ), then the denominator e + (3 ^ is dominated by e, and these 
spurious oscillations cannot be detected by the ESWENO dissipation mechanism. However, if e < 0(5 f 2 ), 
then the weights begin to deviate from their preferred values, which increases dissipation in those stencils 
containing spurious oscillations. 

From these considerations it follows that the lower bound of the constraints on e given by equations 
(62), (72), and (76) should be used (1) to obtain the design order of convergence at the smooth parts of 
the solutions, (2) to effectively damp high-frequency spurious oscillations, and (3) to provide good shock- 
capturing capabilities near strong discontinuities. Therefore, in our numerical experiments, we use 


0(Aa; 3s_4 ) , 

for ip = 0 


0(Ax 3s ~ 3 ) , 

for ip = <p c 

(78) 

0(Ax 3s ~ 2 ) , 

otherwise 



where s — 1 is the degree of the reconstruction polynomials. With the above selection of e, all spurious 
oscillations of amplitude 0(e 1 / 2 ) and higher are suppressed by the ESWENO dissipation mechanism. The 
same conclusions can be drawn for the WENO schemes with the new weights given by equations (58)-(61). 
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VI. Numerical Results 


We now assess the performance of the fourth-, fifth-, and sixth-order ESWENO schemes with the new 
weights and compare them with the conventional WENO counterparts. For all of the numerical experiments 
that are presented, the parameters e of the ESWENO weight functions and the parameters 8, of the additional 
artificial dissipation operator are chosen based on equations (78) and (42) with two minor modifications. 
First, to preserve the scale invariance of the original WENO schemes, the physical grid spacing Aa; in 



Figure 2. Solutions obtained with fourth-order WENO and ESWENO schemes on 101-point grid for linear wave 
equation with initial condition [eq. (81)] at t = 0.04. 


equations (42) and (78) is replaced with A£, where A£ = 1 / J and J is the total number of grid cells. Second, 
as follows from equation (58), the parameter e should be scaled consistently with /3^ r >. In regions where 
the solution is smooth, for the third- and fourth-order WENO and ESWENO schemes, (3' r > approximates 
f' 2 A£ 2 . For the fifth- and sixth-order schemes, (3 ^ approximates the linear combination of the same first- 
order derivative term and a second-order derivative term that is proportional to f" 2 . The same pattern 
persists for higher order schemes as well. In the vicinity of discontinuities, (3 ( r > is of the order 0(/ 2 ). If we 
take into account the above considerations, for all test problems considered the parameter e is set to 

! CA£ 3s_4 , for ip = 0 
CA£ 3s ~ 3 , for ip = (fi c 

CA£ 3s_2 , otherwise (^9) 

C = max (|I/oIUI/o 2 ">--' t II (/o S_1) ) ll) , 

where /o = /[wo(£)] is the initial flux, uq(£) is the initial condition, f^ s ^ is the (s — l)th-order derivative 
of f 0 with respect to £, s — 1 is the degree of the reconstruction polynomials, is a set of points at which 
the solution is discontinuous, and || • || is a norm in which the solution is sought. The scaling factor C in 
equation (79) can easily be evaluated because it depends only on the initial condition for which the locations 
of all discontinuities are known a priori. Note that the parameter e is calculated once and the same value is 
used over the entire time interval of integration; thus the computational cost does not increase. 

In accordance with equation (23), the parameter <5, should be scaled consistently with A;. Because A* is 
a linear combination of the weights with each being of the order 0(1), the scaling factor of S t is set equal to 
one, which leads to 

St = A fP- 2i+1 . (80) 
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Equations (79) and (80) eliminate the ambiguity in determining the parameters e and Si for the ESWENO 
schemes; thus, the equations for e and 8 t are free of tuning parameters. Furthermore, equations (79) and 
(80) are fully consistent with the sufficient conditions [eqs. (42) and (78)] and allow the ESWENO schemes 
to be invariant when the spatial and time variables are scaled by the same factor. Note that the parameter 
e for the conventional WENO schemes of Jiang and Shu is set to 10~ 6 as recommended in. 2 In accordance 
with the recommendations of Borges et ah, 8 the parameter e for the fifth-order WENO scheme with the 
weights given by equations (17)— (19), and (22), which is referred as WENO-Z, is set to 10~ 40 . 

The time derivative for all steady test problems is approximated by using a third-order total variation 
diminishing (TVD) Runge-Kutta method that is developed in reference, 11 while unsteady problems are 
integrated by using a fourth-order low-storage Runge-Kutta method (ref. 12 ). To reduce the fourth-order 
temporal error component and make it consistent with the spatial error of fifth- and sixth-order schemes, the 
time step in global grid refinement studies is reduced by a factor of 2 6 / 4 for each doubling of the number of 
grid points in space. The Courant-Friedrich-Levy (CFL) number has been set to 0.3 and 0.6 for the steady 
and unsteady test problems, respectively. 

VI. A. Scalar Linear Wave Equation 

We begin by verifying that the new class of ESWENO schemes provide the design order of convergence for 
smooth problems, including local extrema. To check this property, we consider equation (1) with a = 1 and 
the following initial condition: 

Uo (x) =e - 300 (*-*e) 2 } (81) 

where x c is 0.5. The computational domain for this test problem is set to 0 < x < 1. Numerical solutions 
are calculated on a sequence of globally refined uniform grids and advanced in time up to t = 1, which 
corresponds to one period in time. 

First, we show that the conventional fourth-order WENO scheme is unstable for this smooth problem; 
however, the corresponding ESWENO is stable, and its solution is in excellent agreement with the exact 
solution, as shown in Figure 2. To make the conventional fourth-order WENO scheme stable, a smoothness 



Figure 3. Loo error norms obtained with the fourth-order linear, WENO, and ESWENO schemes for linear wave 
equation with initial condition given in equation (81). 


indicator that corresponds to the downwind stencil has been modified as follows: 


P R = 


(, 0 L ) k + {j3 c ) k + (f3 R ) k 


l/k 


( 82 ) 
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where k is a constant that is greater than 1. Note that j3 R is a C°° function of its arguments and ap- 
proaches nra x(/3 L , (3 C , f3 R ) as k — > oo. In contrast to equation (82), a modified smoothness indicator 
j3 R = ma x(/3 L , /3 C , /3 R ) that is proposed in reference 6 is a nonsmooth function of (3 L , (3 C and f3 R , which 
may lead to the degeneration of the design order of accuracy and is more prone to spurious oscillations near 
unresolved features and strong discontinuities. For k — > oo, the downwind smoothness indicator that is given 
by equation (82) prevents the corresponding weight w R from being larger than the lesser of the other two 
weight functions; thus, the stencil is biased in the upwind direction near unresolved features. In smooth re- 
gions, all three smoothness indicators are of the same order, and the weights approach their preferred values 
of w L = w c = |, and w R = which provides the design order of accuracy. For all of the numerical 
experiments that are presented herein, the parameter k in equation (82) is equal to 4. 




Figure 4. L 0 0 error norms obtained with the fifth- and sixth-order linear, WENO, WENO-Z, and ESWENO schemes 
for linear wave equation with initial condition given in equation (82). 


Figure 3 shows the L ^ error norms that are obtained with the fourth-order WENO and ESWENO 
schemes, and the corresponding underlying linear schemes. As shown in Figure 3, the fourth-order ESWENO 
scheme is significantly more accurate than the conventional fourth-order WENO scheme. The L ^ error 
that is obtained with the fourth-order ESWENO scheme is slightly higher than that obtained with the 
corresponding linear scheme on coarse meshes and reaches its theoretical limit starting at J = 201 grid 
points. The maximum error occurs at the peak of the Gaussian pulse indicating that the ESWENO scheme 
is design-order accurate at the smooth extrema. In contrast to the ESWENO scheme, the conventional 
WENO scheme demonstrates only a third-order convergence rate, even on the finest mesh with J = 1601, 
and is two to three orders of magnitude less accurate on moderate and fine meshes. 

If the new weights [eqs. (58) and (59)] are used instead of their conventional counterparts, then the design 
order of convergence of the fourth-order central WENO scheme is recovered, as shown in figure 3. Moreover, 
the fourth-order WENO scheme with the new weights provides slightly better accuracy on coarse meshes 
than the corresponding ESWENO scheme. This result is not surprising because the ESWENO scheme has 
the additional dissipation term [eq. (24)], which guarantees the stability. In general, if a WENO scheme with 
the new weights is stable, then it is slightly less dissipative than the corresponding ESWENO counterpart. 
For all of the problems that are considered, however, the results obtained with the conventional WENO 
schemes with the weights given by equations (58) and (59) are practically indistinguishable from those of 
the corresponding ESWENO schemes; therefore these results are not presented hereafter. 

The L x error norms that are obtained with the fifth-order WENO-Z scheme and the fifth- and sixth- 
order WENO and ESWENO schemes for the same test problem are depicted in figure VI. A. Similar to the 
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fourth-order case, the conventional sixth-order central WENO scheme is unstable. This problem is avoided 
by using the same upwinding technique that is outlined earlier. For the sixth-order central WENO scheme, 
the modified smoothness indicator that corresponds to the most downwind stencil is given by 


Both the fifth-order WENO-Z and fifth- and sixth-order ESWENO schemes are equal in accuracy to the 
corresponding underlying linear schemes, as shown in Figure VI. A. Although the fifth-order WENO scheme 
exhibits the design-order convergence rate, its L ^ error norm is nearly an order of magnitude larger than 
those of the corresponding fifth-order WENO-Z and ESWENO schemes. In contrast to that of the sixth-order 
ESWENO scheme, the error that is obtained with the conventional sixth-order WENO scheme shows 
only fifth-order convergence and is quite far from the theoretical limit (represented by the corresponding 
sixth-order central linear scheme). Note that the WENO-Z schemes of fourth- and sixth-order are currently 
unavailable in the published literature and therefore are not presented herein. 

5th-order WENO 


5th-order ESWENO 



Time 

Figure 5. Time histories of rightmost eigenvalue of symmetric part of fifth-order WENO and ESWENO operators for 
Gaussian pulse problem. 

As we have proven in sections 2 and 4.1, for the linear convection equation (1) with piecewise continuous 
initial conditions, the family of ESWENO schemes is stable in the energy norm. The sufficient condition for 
stability is the negative semidefiniteness of —\{D + D T ), where D is the ESWENO discrete operator. This 
property implies that all eigenvalues of the symmetric portion of the ESWENO operator are nonpositive. In 
contradistinction to the ESWENO scheme, the symmetric portion of the WENO operator may have positive 
eigenvalues. (See section 4.1.) Note that the same conclusion can be drawn for the fifth-order WENO-Z 
scheme (ref. 8 ), whose derivative operator is identical to that of the conventional WENO scheme with the 
modified weights [eqs. (17)— (19) , and (22)] which vary in the interval from 0 to 1 near the unresolved features. 
These properties of the WENO, WENO-Z, and ESWENO schemes are shown in figure 5, which gives the 
time histories of the rightmost eigenvalue of the symmetric portion of the operators computed on a 201-point 
grid. The rightmost eigenvalue of the ESWENO operator — ^(D + D T ) is equal to zero up to the order of 
the round-off error, while the symmetric part of the conventional fifth-order WENO and WENO-Z operators 
have positive eigenvalues of order O(10 -1 ) and 0(10), respectively. (See figure 5.) This results from the fact 
that the symmetric portion of these WENO-type operators is not negative semidefinite, if unresolved features 
exist in the computational domain. Note that the presence of the positive eigenvalues does not imply that 
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the fifth-order WENO and WENO-Z schemes are globally unstable because these eigenvalues correspond to 
different grid points at different moments of time. 



Figure 6. L 0 0 error norms obtained with fifth-order linear, WENO, WENO-Z, and ESWENO schemes for linear wave 
equation with initial conditions given in equation (84) at t = 1. 


The superiority of the ESWENO schemes with the new weights compared with the conventional and 
modified WENO counterparts, becomes evident when a more challenging problem is considered. In the 
previous test problem, the initial condition [eq. (81)] has the critical point at which /' = 0 and /" ^ 0 
(i.e. , the number of vanishing derivatives is n v d = 1). As shown in section 5.3, the fifth-order WENO-Z 
scheme that was developed by Borges et al. 8 fails to recover the design order of convergence if n v d > 6 for 
an arbitrary choice of the parameter m in equation (74). We now demonstrate this property for n v d > 3. 
Consider equation (1) with the following initial condition: 

f z 18 - 14 z 16 + 69z 14 - 175z 12 + 259z 10 

u 0 (x) = l -231z 8 + 119z 6 -29z 4 + 1, for \z\ < 1 (84) 

I 0 otherwise 


where z = 5(x — 0.5) and 0 < x < 1. The above function is six times continuously differentiable and has 
three critical points: one at x = 0.5 with n v d = 3 (i.e., /'( 0.5) = /"(0.5) = /'"(0.5) = 0 and /""(0.5) ^ 0) 
and the remaining two at x = 0.3 and 0.7 with n v d = 6. If we compare the L ^ error norms computed 
with the fifth-order linear, WENO, WENO-Z, and ESWENO schemes at t = 1, (see fig. 6), we see that 
the situation changes dramatically as compared with the previous test problem. As expected, the WENO- 
Z scheme fails to deliver fifth order convergence, while the ESWENO scheme is fifth-order accurate and 
provides practically the same error convergence as the underlying linear scheme. This result is not surprising 
because the parameter e for the WENO-Z scheme is set to 10 -40 , as suggested in reference. 8 As a result, 
at a point x = 0.5 + Ax the weights become of order 0(1). Thus, if we use equations [(65), (22), and (74)] 
with e = 0 and take into account that for the polynomial in equation (84), /' = 0( Aa; 3 ), f" = 0(Ax 2 ), 
f" = O(Ax), and f"" = 0(1) at x = 0.5 + Ax. then we have 


w 


X) _ d(0 = O 


\ff"f"~f'f""\Ax 5 
f 2 Ax 2 + (f /' 2 - f'f'")Aa 



(85) 


As follows from the above estimate, the WENO-Z scheme fails to recover fifth order convergence regardless 
of the choice of the parameter m in equation (74). Our numerical results corroborate this conclusion. Figure 
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6 shows that increasing the exponent m in equation (74) results in an even larger deterioration in accuracy. 
In contrast to the WENO-Z scheme, the conventional fifth-order scheme of Jiang and Shu 2 converges at the 
design order and begins to approach the theoretical limit on the finest meshes. The main reason for such 
a behavior is the choice of the parameter e, which is set to 10 -6 . As the grid is refined, f^ r > — > 0, the 
denominator in equation (21) is dominated by e, and the weights approach their preferred values. 

VI. B. The 1-D Euler Equations 

In reference, 1 the third-order ESWENO scheme is proved to be energy-stable for a system of linear hyper- 
bolic equations with periodic boundary conditions. This result can be directly extended to the class of higher 
order ESWENO schemes that are presented in this paper. For nonlinear conservation laws with nonperiodic 
boundary conditions, a similar proof for the high-order ESWENO schemes is not currently available. Never- 
theless, we would like to test how the new schemes perform for the system of the quasi-l-D Euler equations, 
which are given by 





dU | dF . 
dt ' dx _ 

= G 



p 


pu 


pu 

u = 

pu 

E 

, F = 

pu 2 + P 
(E + P)u 

0 

II 

pu z 

_ (E + P)u _ 


P=( 7 -l)(E+*£) , 


where A = A(x) is the cross-sectional area of a quasi-l-D nozzle and 7 = 1.4. As has been shown in 
reference, 1 the ESWENO reconstruction must be implemented in local characteristic fields to guarantee the 
stability of the ESWENO scheme for systems of hyperbolic conservation laws. In the present analysis, both 
the WENO and ESWENO reconstructions are based on the Lax-Friedrichs flux splitting. See reference 1 for 
further details. 


3rd-order linear 5th-order linear 




Figure 7. Loo error norms obtained with the third- and fourth-order and fifth- and sixth-order linear, WENO, and 
ESWENO schemes for the subsonic quasi-l-D nozzle problem. 


For the 1-D Euler equations, the preferred biasing of the stencil in the upwind direction, given by 
equations (82) and (83), is used for the central fourth- and sixth-order WENO and ESWENO schemes to 
suppress spurious oscillations near strong discontinuities. In the ESWENO case, these oscillations can also 
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be eliminated by increasing the coefficient in front of the DiK\D f term in the artificial dissipation operator. 
However, this approach is more dissipative and also reduces the maximum CFL number for which the scheme 
remains stable. For the test problems that are presented in this section, the results obtained with the fifth- 
order WENO-Z scheme are practically indistinguishable from those of the fifth-order ESWENO scheme; 
therefore the WENO-Z results are not presented hereafter. 

To verify that the new ESWENO schemes are design-order accurate for hyperbolic systems, we consider 
the steady-state isentropic flow through a quasi-l-D nozzle with the following cross-sectional area: A(x) = 
1 — 0.8x(l — x), 0 < x < 1, as a test problem. The inflow Mach number is set to 0.5 and the pressure at 
x = 1 is assumed to be equal to that at x = 0. Under these conditions, the flow is fully subsonic, and the 
solution is smooth. Global grid-refinement studies for the third-, fourth-, fifth-, and sixth-order WENO and 
ESWENO schemes are presented in figure 7. 

The Loo error norms that are obtained with the third- and fourth-order ESWENO schemes exhibit the 
design order of convergence. On the finest mesh with J = 801, the fourth-order ESWENO solution is 
approximately three and four orders of magnitude more accurate than those of the third-order ESWENO 
and WENO schemes, respectively. As shown in figure 7, the fifth-order ESWENO scheme provides the 
same accuracy as its underlying linear scheme, thereby exhibiting the perfect error convergence for this test 
problem. Although the conventional fifth-order WENO scheme exhibits the design order of convergence, 
its Loo error norm is larger by a factor of 3 than that obtained with ESWENO counterpart. Similar to 
the fourth-order case, the central sixth-order ESWENO scheme converges at the design-order rate on both 
coarse and fine meshes, thus providing significantly more accurate solutions than both fifth-order schemes. 

Remark 9. For this problem, the conventional central fourth- and sixth-order WENO schemes do not converge 
to a steady-state solution even with the upwinding mechanisms [eqs. (82) and (83)] turned on. The primary 
reason for this tendency is the presence of positive eigenvalues in the spectrum of the symmetric part of the 
conventional WENO operator. 



Figure 8. Comparison of fifth-order WENO and ESWENO schemes for steady transonic flow through quasi-l-D nozzle. 

The next test problem is the steady transonic flow through a quasi-l-D nozzle with the following cross- 
sectional area: 

A{x) = 1.398 + 0.347 tanh(0.8a: — 4), 0 < x < 10. 

The Mach number at x = 0 is 1.5, and the outflow conditions have been chosen so that the shock is located 
at x = 5. Density profiles are calculated using the fifth-order WENO and fifth- and sixth-order ESWENO 
schemes on a 51-point grid and are compared with the exact solution in figure 8. For all schemes, the 
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Density 


captured shock is smeared over two grid cells, and the numerical solutions are essentially nonoscillatory and 
agree quite well with the exact solution. Note that the fifth- and sixth-order ESWENO schemes converge 
to the steady-state solution an order of magnitude faster than the conventional fifth-order WENO scheme; 
this result again indicates the presence of unstable modes that are generated by the positive eigenvalues of 
the WENO dissipation operator. 




Figure 9. Density profiles computed with the fourth-order and sixth-order WENO and ESWENO schemes on a uniform 
grid with 201 points for the Sod problem at t = 0.16. 


The last two problems considered are standard unsteady problems for testing shock-capturing schemes. 
The first problem is the Riemann problem with the initial conditions proposed by Sod: 13 


(p,u,P) = 


( 1 , 0 , 1 ) 

(0.125,0,0.1) 


if - 0.5 < x < 0 
if 0 < x < 0.5 


The numerical solutions that are computed with both fourth- and sixth-order WENO and ESWENO schemes 
are essentially nonoscillatory, as one can see in figure 9. Note, however, that the fourth- and sixth-order 
ESWENO schemes provide better resolution near the contact discontinuity and the shock as compared with 
the conventional WENO schemes. 

We conclude this section by presenting the numerical results that are obtained with WENO and ESWENO 
schemes for the shock entropy-wave interaction problem. The solution of this benchmark problem contains 
both strong discontinuities and smooth structures, and is well suited for testing high-order shock-capturing 
schemes. The governing equations are the time-dependent 1-D Euler equations [eq. (86)] with G = 0, 
subject to the following initial conditions: 


(3.857134, 2.629369, 10.33333) , 
(1 + 0.2 sin 5a;, 0, 1) , 


if — 5 < x < — 4 

if — 4 < x < 5 . 


(87) 


The governing equations are integrated in time up to t = 1.8. The exact solution to this problem is not 
available. Therefore, a numerical solution that is obtained with the conventional fifth-order WENO scheme 
on a uniform grid with J = 4001 grid points is used as a reference solution. 

Numerical solutions that are computed with the fourth-, fifth-, and sixth-order WENO and ESWENO 
schemes on a 301-point grid at t = 1.8 are compared with the “exact” reference solution in figure 10. Both 
the WENO and ESWENO solutions are free of spurious oscillations. The fourth-order WENO scheme is 
the most dissipative among the considered schemes and has the worst resolution in a region just upstream 
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of the moving shock. The fourth-order ESWENO scheme provides better resolution of the high-frequency 
oscillations behind the shock. However, the wave amplitudes are much lower than those that are computed 
with the fifth- and sixth-order schemes. Again, the fifth- and sixth-order ESWENO solutions are more 
accurate than those of the corresponding WENO schemes. Both ESWENO schemes resolve the smooth 
extrema quite well, and the solutions are essentially nonoscillatory near the captured shock. 

VII. Conclusions 

A systematic methodology is developed to construct Energy Stable weighted essentially nonoscillatory 
(ESWENO) finite-difference schemes of arbitrary order, and is used to construct ESWENO schemes based 
on the WENO schemes that are presented in references 2 and . 6 The new ESWENO schemes differ from the 
conventional WENO schemes by the addition of a nonlinear artificial dissipation term of special form. The 
additional term is design-order accurate for smooth solutions, including smooth extrema, and guarantees that 
the ESWENO scheme is stable in the L 2 -energy norm for both continuous and discontinuous solutions of 
hyperbolic systems; for the conventional WENO scheme, an energy estimate is not available. The distinctive 
feature of the new class of ESWENO schemes is that the spectrum of the symmetric part of the ESWENO 
operator is always located in the left half-plane, while the symmetric part of the WENO operator has positive 
eigenvalues; thus, the conventional WENO schemes may become locally unstable. 

New weight functions are also developed, guided by newly derived sufficient conditions guaranteeing 
design order accuracy for arbitrary WENO schemes; sufficient conditions that are valid regardless of the 
number of vanishing solution derivatives. A truncation error analysis is used to optimize the weight functions, 
so that the new ESWENO schemes satisfy the sufficient conditions for accuracy and provide essentially 
nonoscillatory solutions for problems with strong discontinuities. Numerical experiments confirm that the 
high-order ESWENO schemes with the new weights are significantly more accurate than conventional WENO 
schemes of the same order, and demonstrate excellent shock-capturing capabilities. 

VIII. Appendix 

VIII. A. Existence of Symmetric Decomposition 

The following lemma establishes the existence of the decomposition D sym = [D\ *] A,; [D\ l ] . The lemma 

provides an algorithm to decompose a matrix of arbitrary bandwidth 2 s + 1 into the sum of two matrices: 
one that can be expressed as one term of the desired symmetric form, and one of bandwidth 2 (s — 1 ) + 1 . 

Lemma 1. Define E s to be an iV-dimensional, symmetric, bidiagonal matrix with the subdiagonal elements 
e s +i.i, e s _|_ 2 , 2 , ••• , ejv-i,Ar-s-i, Cn,n-s- The parameter s is the offset from the main diagonal. Further, 
define E s _i to be an TV-dimensional, symmetric, banded matrix of bandwidth 2(s — 1) + 1 with elements that 
are functions of the matrix E s . The matrices E s and B s _ \ satisfy the following relationship: 

E s = (— 1) S [D 1 S ]A S [D 1 S ] T + E S _ 1 . 

with A s = diag[e s +i t i, e s + 2 , 2 , • • • , ejv-i.jv-s-i, £n,n-s, 0 s t ] and 0 S is a zero vector of dimension s. 

Proof: A proof by induction is not included herein. One step of the algorithm can be verified immediately 
by inspection. 

A simple example illustrates Lemma 1. Define the 5x5 matrices E 2 , D 15 , and B\ such that E 2 = 
(— l) 2 [I?i 5 "]A 2 [Di 5 2 ] + Bi and 
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/ 
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a 

—2a 
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\ 
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b 
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0 




—2a 

4 a + b 

— 2(a + b) 
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A 2 = 


0 

0 

c 
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0 


; B 1 = 


0 

—2 (a + b) 

a + Ab + c 

— 2(6 + c) 
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0 

0 




0 

0 

—2(6 + c) 

6 + 4c 

-2c 



V 

0 

0 

0 

0 

0 

J 



0 

0 

0 

-2c 

c 

) 


Remark 10. Lemma 1 can be used recursively, beginning with the outermost s diagonal and working inward, 
to establish that D sym = ■ Each step generates a new symmetric matrix with a bandwidth 

that is equal to 2 less than the previous value until only a diagonal matrix remains. Thus, the algorithm 
terminates after s + 1 steps with the desired decomposition. 

Remark 11. This recursive algorithm works for nonperiodic matrices in multiple spatial dimensions. 


VIII. B. ESWENO Schemes of Third- and Fourth-Order 

Based on equation (13), / J+ i defined for a third- or fourth-order WENO scheme is 



= w 


L fL 
j+l/lJ 3 + 1/2 


W? 


c 


j+l/2Jj+l/2 


+ W 


R 


R 


j+l/2Jj+l/2 


( 88 ) 


where w L , w c , and w R are weight functions that are assigned to each respective stencil. The second-order 
linear fluxes: /< / 2 , defined in equation (88) for r = { L , C, R }, are 


( f L (u j+ 1 /2 ) \ 

f C i U 3 + 1/2) j 

V f R (u j+ 1 / 2 ) ) 



-1 

0 

0 


3 0 0 

1 1 0 
0 3-1 


( f(uj- 1 ) \ 
f(Uj) 
f(u j+ 1 ) 

V /( u i+ 2) / 


(89) 


The terms w L , w c 
with 


and w R are the nonlinear weight functions that are given by equations (58) and (59) 


+3 = 


(~fj- 1 + 3 fj - 3/j+i + /j+ 2 ) 2 , for <p ^ 0 
(fj- 1 - 2/j + fj+ 1 ) 2 , for <p = 0 


(See ref. 1 for further details). These weight functions have preferred values that are derived from underlying 
linear schemes, as well as solution-dependent components. The preferred values are given by the formulae 


d R = | — ip ; dP = I 


d R = 




(90) 


where ip is a parameter. The convergence rate of equation (13), with the preferred weight values that are 
given in equation (88), is equal to 3 for all values of the parameter ip except p c = for which the convergence 
rate is 4. 

Combining equation (89) with equation (88) and substituting the resulting WENO flux / J+ 1 into equation 
(13) produces a stencil for the jth grid point of the form 



( 0 

W j + 1/2 

3%'+ 1/2 

0 

° ^ 
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0 

<+1/2 

<+1/2 
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1 
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0 

3 <+i/2 

-<+1/2 

2Ax 

W f-l/2 

-3uP- 1/2 

0 
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0 

- w ° 3 - 1/2 

~ w f - 1/2 

0 

0 


l 0 

0 

~ 3w f-l/2 

<-1/2 

0 ) 


with k on the interval j — 2 le k < j + 2. 


/ f(Uj- 2 ) \ 

f(uj- 1 ) 
f(Uj ) 
f{Uj+l) 

V f(u j+ 2 ) / 


( 91 ) 
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Simplifying D 3 using the expression 


w c = 1 - w L - w R 


yields expressions for a single row of the D 3 matrix of the following form 

/ : 

0 


D 3 = 


j ’ k 2Ax 


W f~ 1/2 

„L 


~ 2w f-l/2 - W^ + l/2 + <-1/2 - 1 

W j- 1/2 + 2u, j+l/2 ~ ^ W j~i/2 — W j+l/2 
“<H/2 + <-1/2 + 2 <+l/2 + 1 
-<+1/2 


V 


( 92 ) 


) 


D 3 

^ skew 


The derivative matrix D 3 is now decomposed into symmetric and skew-symmetric parts 

r)3 _ r)3 , rj3 

J ^ skew ' ^ sym 

As with the fifth-order case, the skew-symmetric component and the norm of D 3 take the forms 

= P-'Qs ; Qs + Ql = 0 ; P = Axl . 
while the matrix D 3 ym is expressed as 

Km = P ~ l ( D ' A i [^1 2 ] T + D l l A ? [Vl 1 ? + A 0 [Dl°] T ) 

where A 3 is a diagonal matrix with the expressions for the jth element defined by 

i r- 

„R 


( A <b 


1 

4 L 


<+ 3/2 W j+ 1/2 


<*?>« = I L 


-<+ 3/2 + <+ 1/2 - <- 1/2 + <- 1/2 


( A 0 < = 2 I -<- 1/2 - <- 1/2 - <- 


J — 1/2 j — 1/2 


J — 1/2 _ 


+<+ 1/2 +<+ 1/2 +<+ 1/2 


= 0 . 


(93) 

(94) 

(95) 

(96) 


If we define (Af) , * = 1, 2, to be smoothly positive, then 


« = JkA a ?)L+*? - («> 


2 L V v 2 

and the additional artificial dissipation operator is given by 

2 


Did = p-^d^DIDi 1 } 


The resulting energy-stable scheme is obtained by adding the additional dissipation term to the corresponding 
conventional WENO scheme as 

D 3 = D 3 + D 3 ad . 
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VIII. C. ESWENO Schemes of Seventh- and Eighth-order 

The fj + 1, defined for the seventh- and eighth-order WENO schemes, is 


fj+i = w 


LL 


LL 


j+l/2Jj+l/2 


W R 


j+l/2fj+l/2 
RR fRR 


+ W 


R 


R 


j+l/2Jj+l/2 


+ W 


■ w 


c 


c 


j+l/2Jj+l/2 


j+l/2Jj+l/2’ 


(97) 


where w LL , w L , w G , w R , and w RR are weight functions that are assigned to each respective stencil, and 
are given by equations (58) and (59) with 


T7 = 


_ J (/,•_ 3 - lfj -2 + 21/, -i - 35 /, + 35/, +1 - 21 /, +2 + 7/, +3 - /, +4 ) 2 , for p ± 0 


, (~fj- 3 + 6/,_2 - 15/, -i + 20 fj - 15/, +1 + 6/, -|_2 - /, +3 ) 2 , for <p = 0 

The fourth-order linear fluxes /j+i/2 for r = {LL, L, C, R, RR}, which are used in equation (97) are defined 


as 


/ -3 13 -23 25 0 
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( f(uj- 3 ) \ 
f(uj- 2 ) 
f(uj- 1 ) 

f(Uj) 
f(.Uj+ 1 ) 

/K+2) 

f(u j+ 3 ) 

V /K+4) / 

The derivative matrix I?' is now decomposed into symmetric and skew-symmetric parts as 


( f LL (u j+ 1 / 2 ) \ 
f L (Uj+ 1 / 2 ) 
f c (uj+ 1 / 2 ) 

f R (u j+1/2 ) 

V /^(M, + l/2) ) 


1 0 


0 25 -23 13 -3 / 


(98) 


ir = d‘ 


D[ 


7 skew 1 sym 

As with the fifth-order case, the skew-synnnetric component of D 7 takes the form 

D 7 skew = P- l Q 7 ; Qr + Qf = 0 ; P = Axl . 


The scheme converges with seventh-order accuracy if the weights are equal to their preferred values 


d LL = ±-p,d L = ±§ - 8 p,d G = § ,d R = ± + 8 p,d RR = p 

and with eighth-order accuracy for the specific value ip c = 

The D 7 sym can be simplified into the form 


D 7 

sym 


= p 


(d^ A l [D^f + Di 3 A 7 [Dd 3 ] 7 ' + Hi 2 A l [D^f 
+ dSaUd, 1 ] 7 + D 1 0 A 7 [D 1 0 f) 


where A ' is a diagonal matrix with expressions for the jth element defined by 




m LL - in RR 
W j+ 7/2 W j+ 1/2 


(^3 — 24 I h™j+ 5/2 9u b+7/2 + 


,'+5/2 


j+7/2 a ','+5/2 


-«f +1 /2 + 9 <-l/2 - *1/2 


(99) 


( 100 ) 


( 101 ) 

( 102 ) 
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Ni)jj 24 I +2lUj + 3/2 llWj +5 / 2 + 9Wj +7 / 2 


+2w j+3/2 ^ W j+ 5/2 


— 10 


C 

3 + 1/2 


- tO 


c 

3+3/2 


( A l)i,7 = 


+ 2 <— 1/ 2 - 2w f+l/2 
-9<4/2 + H<-l/2 - 2 <+l/2 


+6w £+ 1/2 ~ 13u7 i+3/2 + 1 ° W f+ i 5/2 - HV 
+ 3w; i+ 1/2 — ^ W i+3/2 "I" W i+5/2 
+^7-1/2 - «>i+3/2 


— W 


R 

i- 3/2 


+ 4 w ^_ 1(/2 - 3to 


i+1/2 


-T^i- 5/2 1U %-3/2 + 10U7 i-l/2 0U? i+l/2 


( A o)j,j - 3 [ 


—wt% , 0 — wf n , 0 — tor i /o — to;* -i /o — w. 


+w. 


LL 
'i- 1/2 
LL 


7+1/2 


u i-l/2 

„ L 

U i+l/2 


, f 

U i- 1/2 
' ^7+1/2 


„R 
i- 1/2 


iti? 

7 — 1/2 

it it 


' 1/2 + W 7+l/2 


= 0 


If we Define (Aj) • , t = 1,4, to be smoothly positive as 


( A 7 ? )j,j - p[\/ ( A DL +<5 7 2 - ( A i )/./] > 


the additional artificial dissipation operator is given by 

4 


i=l 


( 103 ) 


(104) 


(105) 


and the resulting energy-stable scheme is obtained by adding the additional dissipation term to the corre- 
sponding conventional WENO scheme as 


D 7 = D 7 


D 


ad 
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